WCM Biostatistics Computing Club: Plotly Tutorial

Devil’s Path Hiking Data

Today we will be working with data from the Devil’s Path, a hiking trail in the Catskill Mountains of New York. It is supposed to be the toughest hike in the Catskills mountains! The data that I accessed is from CalTopo (https://caltopo.com/), a mapping and trip planning tool for the backcountry. If you are a member of CalTopo you can get a .csv of any hiking trail’s latitude, longitude, distance, and elevation. I have provided the data from CalTopo in the github repository.

devils.path <- read.csv('Devils_Path.csv') 
knitr::kable(head(devils.path))

Lat	Lng	Distance..meters.	Distance..feet.	Distance..miles.	Elevation..meters.	Elevation..feet.	Slope..degrees.	Aspect..degrees.	Landcover	Canopy..percent.
42.13400	-74.10404	0	0	0.0000000	626	2054	9	40	Forest	79
42.13383	-74.10510	90	294	0.0556631	632	2073	6	339	Forest	85
42.13324	-74.10629	207	681	0.1288833	628	2060	4	310	Forest	89
42.13232	-74.10648	310	1018	0.1928586	630	2067	3	308	Forest	86
42.13137	-74.10583	430	1410	0.2669656	633	2077	8	192	Forest	88
42.13041	-74.10542	542	1779	0.3368528	640	2100	11	304	Forest	90

What is Plotly?

• Plotly is an interactive, browser-based charting library built on the open source javascript graphing library, plotly.js. It works entirely locally, through the HTML widgets framework.

• Plotly has interfaces with R, Python, Matlab and Julia

• Plotly plays well with Rmarkdwon! If you are using Rmarkdown with HTML output, printing a plotly object in a code chunk will result in an interactive HTML graph. When using rmarkdown with non-HTML output, printing a plotly object will result in a png screenshot of the graph.

library(plotly)

Plotly: Our Example

Here is the plot that we will be working towards making today!

The Elevation Plot

Here is a simple plot of the distance on the trail versus the elevation, using the plot_ly function from the plotly package.

plot_ly(data = devils.path,
        x = ~Distance..miles.,
        y = ~Elevation..feet.)

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

If you are annoyed by the warning messages you can specify the type and the mode in the plot_ly command (but honestly the command does a pretty good job of figuring out what you are after!)

plot_ly(data = devils.path,
        x = ~Distance..miles.,
        y = ~Elevation..feet.,
        type = 'scatter', 
        mode = 'markers') 

You can change the plot to a line plot by specifying mode = ‘lines’ instead of mode = ‘markers’ (the default). You can also change the color of the line.

plot_ly(data = devils.path,
        x = ~Distance..miles.,
        y = ~Elevation..feet.,
        type = 'scatter', 
        mode = 'lines',
        color = I('purple')) 

Changing the axis titles in the plot is a little cumbersome… the title of the plot itself is easier.

x <- list(
  title = "Distance (in miles)"
)

y <- list(
  title = "Elevation (in feet)"
)

plot_ly(data = devils.path,
        x = ~Distance..miles.,
        y = ~Elevation..feet.,
        type = 'scatter', 
        mode = 'lines',
        color = I('purple')) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Devil's Path")

You can also change what the hover on the plot says using ‘text’ and ‘hoverinfo’.

x <- list(
  title = "Distance (in miles)"
)

y <- list(
  title = "Elevation (in feet)"
)

plot_ly(data = devils.path,
        x = ~Distance..miles.,
        y = ~Elevation..feet.,
        type = 'scatter', 
        mode = 'lines',
        color = I('purple'),
        hoverinfo = 'text',
        text = paste(
          round(devils.path$Distance..miles., 1), 
          "miles, ",                    
          round(devils.path$Elevation..feet., 0),
          "feet elevation")) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Devil's Path")

It’s also a big deal for mountain’s in the Catskills to be over 3500 feet, so we will draw a horizontal line at 3500 feet on our plot. We will also fill in the plot and take it all the way down to 0 elevation.

x <- list(
  title = "Distance (in miles)"
)

y <- list(
  title = "Elevation (in feet)"
)

plot_ly(data = devils.path,
        x = ~Distance..miles.,
        y = ~Elevation..feet.,
        type = 'scatter', 
        mode = 'lines',
        color = I('purple'),
        hoverinfo = 'text',
        text = paste(
          round(devils.path$Distance..miles., 1), 
          "miles, ",                    
          round(devils.path$Elevation..feet., 0),
          "feet elevation"),
        fill = 'tozeroy')  %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Devil's Path",
         shapes=list(type='line', 
                     y0= 3500, 
                     y1= 3500, 
                     x0=min(devils.path$Distance..miles.),
                     x1=max(devils.path$Distance..miles.),
                     line=list(dash='dot', width=1)))

A Digression into Exploring Some of the Other Variables…

What if we want to take a look at the elevation versus the percent canopy cover? And what exactly is canopy cover? Canopy cover is the percentage of total ground area that is covered by the vertical projection of tree crowns. We would expect the canopy cover to decrease as the elevation increases.

plot_ly(data = devils.path,
        x = ~Elevation..feet.,
        y = ~Canopy..percent.)

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

That is the exact trend that we see! Let’s put a title and label the axes of the plot.

x <- list(
  title = "Elevation (in feet)"
)

y <- list(
  title = "Canopy Percentage"
)

plot_ly(data = devils.path,
        x = ~Elevation..feet.,
        y = ~Canopy..percent.) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Devil's Path \n Canopy Percentage versus Elevation")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

And let’s add color by the type of landcover (most of this is forest so this is not visually very interesting – but it is good to know how to do it).

x <- list(
  title = "Elevation (in feet)"
)

y <- list(
  title = "Canopy Percentage"
)

plot_ly(data = devils.path,
        x = ~Elevation..feet.,
        y = ~Canopy..percent.,
        marker = list(
          color = factor(devils.path$Landcover,labels=c("brown",
                                                        "green",
                                                        "blue"))  
          )) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Devil's Path \n Canopy Percentage versus Elevation")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

x <- list(
  title = "Elevation (in feet)"
)

y <- list(
  title = "Canopy Percentage"
)

plot_ly(data = devils.path,
        x = ~Elevation..feet.,
        y = ~Canopy..percent.,
        hoverinfo = 'text',
        text = paste(
          'Canopy Percentage: ', 
          devils.path$Canopy..percent., 
          '%',
          '\n', 
          'Elevation: ',                    
          devils.path$Elevation..feet.,
          ' feet \n',
          'Landcover: ',
          devils.path$Landcover),
        marker = list(
          color = factor(devils.path$Landcover,labels=c("brown",
                                                        "green",
                                                        "blue"))  
          )) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Devil's Path \n Canopy Percentage versus Elevation")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

A taller mountain may be more interesting to look at canopy percentage versus elevation. Mt. Bierstadt is a 14er in Colorado (elevation of above 14,000 feet!). Let’s look at the data…

bierstadt <- read.csv('Bierstadt.csv') 
knitr::kable(head(bierstadt))

Lat	Lng	Distance..meters.	Distance..feet.	Distance..miles.	Elevation..meters.	Elevation..feet.	Slope..degrees.	Aspect..degrees.	Landcover	Canopy..percent.
39.59713	-105.7100	0	0	0.0000000	3539	11611	4	34	Shrub	0
39.59697	-105.7099	18	58	0.0110687	3540	11614	4	61	Shrub	0
39.59681	-105.7099	36	117	0.0221374	3541	11617	4	68	Shrub	0
39.59668	-105.7100	51	167	0.0315679	3541	11617	5	68	Shrub	0
39.59658	-105.7101	69	225	0.0426366	3543	11624	6	57	Shrub	0
39.59647	-105.7102	85	277	0.0525474	3544	11627	7	50	Shrub	0

And plot the canopy percentage versus elevation.

x <- list(
  title = "Elevation (in feet)"
)

y <- list(
  title = "Canopy Percentage"
)

plot_ly(data = bierstadt,
        x = ~Elevation..feet.,
        y = ~Canopy..percent.,
        hoverinfo = 'text',
        text = paste(
          'Canopy Percentage: ', 
          bierstadt$Canopy..percent., 
          '%',
          '\n', 
          'Elevation: ',                    
          bierstadt$Elevation..feet.,
          ' feet \n',
          'Landcover: ',
          bierstadt$Landcover),
        marker = list(
          color = factor(bierstadt$Landcover,labels=c("brown",
                                                      "orange",
                                                      "green",
                                                      "black",
                                                      "blue"))  
          )) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "Mt. Bierstadt \n Canopy Percentage versus Elevation")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

We can jitter the y-values a little to make the plot more attractive.

x <- list(
  title = "Elevation (in feet)"
)

y <- list(
  title = "Canopy Percentage"
)

plot_ly(data = bierstadt,
        x = ~Elevation..feet.,
        y = ~jitter(Canopy..percent., factor = 10),
        hoverinfo = 'text',
        text = paste(
          'Canopy Percentage: ', 
          bierstadt$Canopy..percent., 
          '%',
          '\n', 
          'Elevation: ',                    
          bierstadt$Elevation..feet.,
          ' feet \n',
          'Landcover: ',
          bierstadt$Landcover),
        marker = list(
          color = factor(bierstadt$Landcover,labels=c("brown",
                                                      "orange",
                                                      "green",
                                                      "black",
                                                      "blue"))  
          )) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "Mt. Bierstadt \n Canopy Percentage versus Elevation")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

We can make a plot to let us know what to anticipate on the trail!

x <- list(
  title = "Distance (in miles)"
)

y <- list(
  title = "Elevation (in feet)"
)

plot_ly(data = bierstadt,
        x = ~Distance..miles.,
        y = ~Elevation..feet.,
        hoverinfo = 'text',
        text = paste(
          'Distnace: ', 
          bierstadt$Distance..miles., 
          ' miles',
          '\n', 
          'Elevation: ',                    
          bierstadt$Elevation..feet.,
          ' feet \n',
          'Landcover: ',
          bierstadt$Landcover),
        marker = list(
          color = factor(bierstadt$Landcover,labels=c("brown",
                                                      "orange",
                                                      "green",
                                                      "black",
                                                      "blue"))  
          )) %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "Mt. Bierstadt")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

The Map Plot

We can start by ploting the latitude and longitude of the path.

plot_ly(data = devils.path,
        y = ~Lat,
        x = ~Lng) 

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

We can even plot this in 3d with plotly (the latitude, longitude and elevation)

plot_ly(data = devils.path,
        y = ~Lat,
        x = ~Lng,
        z = ~Elevation..feet.,
        mode = 'lines') 

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter3d

It would look a lot better if we had this on a map! One way to do this is using the ‘plot_mapbox’ function. To use this function, you first need to make an account on mapbox (https://www.mapbox.com). Once you make an account, click on the dropdown menu in the upper right hand corner, go to ‘Account’ and then down to ‘Access tokens’. Paste your access token in here:

Sys.setenv('MAPBOX_TOKEN' = 'your access token') 

Now we can plot the latitude and longitude of the trail on a map from mapbox.

plot_mapbox(data = devils.path,
            lat = ~Lat,
            lon = ~Lng) %>%
  layout(mapbox = list(zoom = 9,
                       center = list(lat = median(devils.path$Lat),
                                     lon = median(devils.path$Lng))))

## No scattermapbox mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

And we can change the style of the map to ‘outdoors’.

plot_mapbox(data = devils.path,
            lat = ~Lat,
            lon = ~Lng) %>%
  layout(mapbox = list(zoom = 9,
                       style = "outdoors",
                       center = list(lat = median(devils.path$Lat),
                                     lon = median(devils.path$Lng))))

## No scattermapbox mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

Or we can chance the style of the map to ‘satellite-streets’.

plot_mapbox(data = devils.path,
            lat = ~Lat,
            lon = ~Lng) %>%
  layout(mapbox = list(zoom = 9,
                       style = "satellite-streets",
                       center = list(lat = median(devils.path$Lat),
                                     lon = median(devils.path$Lng))))

## No scattermapbox mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

##Putting Everything Together with Cross Talk

Now we are ready to plot the elevation plot and the map plot together. We will use the ‘crosstalk’ R package. Crosstalk is an add-on to the htmlwidgets package. It extends htmlwidgets with a set of classes, functions, and conventions for implementing cross-widget interactions.

We will start with the function ‘SharedData’. The primary use for SharedData is to be passed to crosstalk-compatible widgets in place of a data frame. Each SharedData$new(…) call makes a new “group” of widgets that link to each other, but not to widgets in other groups.

library(crosstalk)

devils.path.sd <- SharedData$new(devils.path)


x <- list(
  title = "Distance (in miles)"
)

y <- list(
  title = "Elevation (in feet)"
)


bscols(
  plot_mapbox(data = devils.path.sd,
                lat = ~Lat,
                lon = ~Lng)  %>%
     layout(mapbox = list(zoom = 9,
                          style = "outdoors",
                          center = list(lat = median(devils.path$Lat),
                                        lon = median(devils.path$Lng)))) %>% 
    highlight(dynamic = TRUE, persistent = TRUE),
  plot_ly(devils.path.sd,     
          x = ~Distance..miles.,
          y = ~Elevation..feet.,
          type = 'scatter',
          hoverinfo = 'text',
          text = paste(round(devils.path$Distance..miles., 1), 
                       "miles, ",                   
                       round(devils.path$Elevation..feet., 0),
                       "feet elevation"))   %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Devil's Path",
         shapes=list(type='line', 
                     y0= 3500, 
                     y1= 3500, 
                     x0=min(devils.path$Distance..miles.),
                     x1=max(devils.path$Distance..miles.),
                     line=list(dash='dot', width=1))) %>% 
    highlight("plotly_selected", persistent = TRUE)
)

The Big Bushwack

big.bushwack <- read.csv('Big_Bushwack.csv') 
knitr::kable(head(devils.path))

Lat	Lng	Distance..meters.	Distance..feet.	Distance..miles.	Elevation..meters.	Elevation..feet.	Slope..degrees.	Aspect..degrees.	Landcover	Canopy..percent.
42.13400	-74.10404	0	0	0.0000000	626	2054	9	40	Forest	79
42.13383	-74.10510	90	294	0.0556631	632	2073	6	339	Forest	85
42.13324	-74.10629	207	681	0.1288833	628	2060	4	310	Forest	89
42.13232	-74.10648	310	1018	0.1928586	630	2067	3	308	Forest	86
42.13137	-74.10583	430	1410	0.2669656	633	2077	8	192	Forest	88
42.13041	-74.10542	542	1779	0.3368528	640	2100	11	304	Forest	90

big.bushwack.sd <- SharedData$new(big.bushwack)


x <- list(
  title = "Distance (in miles)"
)

y <- list(
  title = "Elevation (in feet)"
)


bscols(
  plot_mapbox(data = big.bushwack.sd,
                lat = ~Lat,
                lon = ~Lng)  %>%
     layout(mapbox = list(zoom = 9,
                          style = "outdoors",
                          center = list(lat = median(big.bushwack$Lat),
                                        lon = median(big.bushwack$Lng)))) %>% 
    highlight(dynamic = TRUE, persistent = TRUE),
  plot_ly(big.bushwack.sd,     
          x = ~Distance..miles.,
          y = ~Elevation..feet.,
          type = 'scatter',
          hoverinfo = 'text',
          text = paste(round(big.bushwack$Distance..miles., 1), 
                       "miles, ",                   
                       round(big.bushwack$Elevation..feet., 0),
                       "feet elevation"))   %>%
  layout(xaxis = x, 
         yaxis = y,
         title = "The Big Bushwack",
         shapes=list(type='line', 
                     y0= 3500, 
                     y1= 3500, 
                     x0=min(big.bushwack$Distance..miles.),
                     x1=max(big.bushwack$Distance..miles.),
                     line=list(dash='dot', width=1))) %>% 
    highlight("plotly_selected", persistent = TRUE)
)

We can also make a 3d plot of the Big Bushwack.

plot_ly(data = big.bushwack,
        y = ~Lat,
        x = ~Lng,
        z = ~Elevation..feet.,
        mode = 'lines') 

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter3d

Plotly and ggplot

big.bushwack$trail <- 'Bushwack'
devils.path$trail <- 'Devil'
trail.data <- rbind(devils.path, big.bushwack)

library(ggplot2)

g <- ggplot(trail.data, 
       aes(x = Distance..miles.,
           y = Elevation..feet.,
           color = trail)) + 
  ylab('Elevation (in feet)') + 
  xlab('Distance (in miles)') + 
  ggtitle("The Devil's Path and The Big Bushwack") + 
  geom_line() + 
  theme(legend.title = element_blank())

g

ggplotly(g)

Your Turn!

Checkout the ‘Slide.csv’ or the ‘Bierstadt.csv’ files and make your own interactive plots of Slide Mountain (Catskills, NY) or Mount Bierstadt (Rocky Mountains, CO). Or ask me to download your favorite trail from Caltopo.